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Abstract 

Image registration is one important task in many image processing applications. It aims 
to align two or more images so that useful information can be extracted through comparison, 
combination or superposition. This is achieved by constructing an optimal transformation 
which ensures that the template image becomes similar to a given reference image. Although 
many models exist, designing a model capable of modelling large and smooth deformation 
field continues to pose a challenge. This paper proposes a novel variational model for im¬ 
age registration using the Gaussian curvature as a regulariser. The model is motivated by 
the surface restoration work in geometric processing [Elsey and Esedoglu, Multiscale Model. 
Simul., (2009), pp. 1549-1573]. An effective numerical solver is provided for the model using 
an augmented Lagrangian method. Numerical experiments can show that the new model 
outperforms three competing models based on, respectively, a linear curvature [Fischer and 
Modersitzki, J. Math. Imaging Vis., (2003), pp. 81-85], the mean curvature [Chumchob, 
Chen and Brito, Multiscale Model. Simul., (2011), pp. 89-128] and the diffeomorphic demon 
model [Vercauteren at al., Neuroimage, (2009), pp. 61-72] in terms of robustness and accu¬ 
racy. 

Key words: Image registration, Non-parametric image registration, Regularisation, Gaus¬ 
sian curvature, surface mapping. 


1 Introduction 

Image registration along with image segmentation are two of the most important tasks in imaging 
sciences problems. Here we focus on image registration. Much research work has been done; for 
an extensive overview of registration techniques see jMUHlII]. The methods can be classified 
into parametric and non-parametric image registration based on the geometric transformation. 
For the first category, the transformation is governed by a finite set of image features or by 
expanding a transformation in terms of basis functions. The second category (which is our main 
concern in this paper) of non-parametric image registration methods is not restricted to a certain 
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parametrisable set. The problem is modelled as a functional minimisation via the calculus of 
variations. Given two images, the reference R and template T, the functional consists of a 
distance measure V(T,R,u) and a regularisation term S(u) where u = (ui(x),U 2 (x)) is the 
sought displacement vector at pixel s £ 11 C R 2 . The term S(u) removes the ill-posedness of 
the minimisation problem. We use the squared L 2 norm of the distance measure to quantify the 
differences between T and R as follows 

V{T, R,u) = l I (T(x + u(x)) - R(x)) 2 dn. (1) 

2 J n 

The distance measure in equation m is the sum of the squared difference (SSD) which is com¬ 
monly used and optimal for mono-modal image registration with Gaussian noise. For multi-modal 
image registration where T, R cannot be compared directly, other distance measures must be used 
[8]. Generally, the regularisation terms are inspired by physical processes such as elasticity, dif¬ 
fusion and motion curvature. As such, elastic image registration was introduced in T] which 
assumed that objects are deformed as a rubber band. 

In previous works, higher order regularisation models HEI were found to be the most ro¬ 
bust while the diffeomorphic demon model [22: offers the most physical transform in terms of 
(nearly) bijective mapping. Diffusion and total variation regularisation models based on first 
order derivatives are less complicated to implement but are at a disadvantage compared to higher 
order regularisation models based on second order derivatives due to two reasons. First, the for¬ 
mer methods penalise rigid displacement. They cannot properly deal with transformations with 
translation and rotation. Second, low order regularisation is less effective than high order one in 
producing smooth transformations which are important in several applications including medical 
imaging. The work of HEBE] proposed a high order regularisation model named as a linear 
curvature, which is an approximation of the surface (mean) curvature and the model is invariant 
to affine registration. This work was later refined in unsung where a related approximation 
to the sum of the squares of the principal curvatures was suggested and higher order boundary 
conditions were recommended. Without any approximation to the mean curvature, the works in 
Eiia developed useful numerical algorithms for models based on the nonlinear mean curvature 
regularisation and observed advantages over the linear curvature models for image registration; 
however the effect of mesh folding (bijective maps) was not considered. The diffeomorphic de¬ 
mon model |23j is widely used due to its property of bijective maps; however the bijection is not 
precisely imposed. Another useful idea of enforcing bijection, beyond the framework we consider, 
is via minimising the Beltrami coefficient which measures the distortion of the quasi-conformal 
map of registration transforms [T2] . 

In this paper we propose a high order registration model based on Gaussian curvature and 
hope to achieve large and smooth deformation without mesh folding. Although the Gaussian 
curvature is closely related to the mean curvature, it turns out our new model based on the 
Gaussian curvature is much better. The motivation of the proposed model comes from two 
factors. Firstly, we are inspired by the work of Elsey and Esedoglu [4] in geometry processing 
where Gaussian curvature of the image surface is used in a variational formulation. The authors 
proposed the Gaussian curvature as a natural analogue of the total variation of Rudin, Osher 
and Fatemi (ROF) [19] model in geometry processing. Aiming to generalise the ROF model 
to surface fairing where the convex shapes in 3D have similar interpretation to the monotone 
functions in ID problems for the ROF model, they showed that, based on the Gauss Bonnet 
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theorem, the complete analogy of the total variation regularisation for surface fairing is the energy 
functional of the Gaussian curvature. A very important fact pointed out in [4] stated that the 
mean curvature of the surface is not a suitable choice for surface fairing because the model is not 
effective for preserving important features such as creases and corners on the surface (although 
the model is still effective for removing noise). Their claims are also supported by the work of [l~3j 
where the authors illustrated several advantages of Gaussian curvature over mean curvature and 
total variation in removing noise in 2D images. First, Gaussian curvature preserves important 
structures such as edges and corners. Second, only Gaussian curvature can preserve structures 
with low gradient. Third, the model is effective in removing noise on small scale features. Thus, 
we believe that Gaussian curvature is a more natural physical quantity of the surface than mean 
curvature. Here we investigate the potential of using Gaussian curvature to construct a high 
order regularisation model for non-parametric image registration of mono-modal images. 

The outline of this paper is as follows. In §2 we review the existing models for non-parametric 
image registration with focus on the demon, linear and mean curvature models. In §3 we introduce 
the mathematical background of Gaussian curvature for surfaces. In §4 we introduce a Gaussian 
curvature model and a numerical solver to solve the resulting Euler-Lagrange equations. We show 
in §5 some numerical tests including comparisons. Finally, we discuss the parameters’ selection 
issue for our model in §6 and present our conclusions in §7. 

2 Review of Non-parametric Image Registration 

In image registration, given two images T and R (which are assumed to be compactly supported 
and bounded operators T, R : C R d — > R + ), the task is to transform T to match R as closely 
as possible. Although we consider d = 2 throughout this paper, with some extra modifications, 
this work can be extended to d = 3. In non-parametric image registration, the transformation is 
denoted by where ip is a vector valued function 

<p(x) : R 2 —» R 2 

where x = (x,y). To separate the overall mapping from the displacement field, we will define 

<p{x) = x + u(x) 

where u(x) is the displacement field. Thus, finding u(x) is equivalent to finding <p. A non- 
parametric image registration model takes the form 

min^ r (u(®)) = T>(T, R, u(x)) + jS(u(x)) (2) 

u(x) 

where the distance measure T> is given as in m and the choice of regulariser S(u(x)) differentiates 
different models. Here u(x) is searched over a set U of admissible functions that minimise J 1 in 
©■ Usually, the set U is a linear subspace of a Hilbert space with Euclidean scalar product. 

The force term /(it), to be used by all models, is the gradient of ©) with respect to the 
displacement field u{x) 

/(«) = ( fi(u ), / 2 (m)) T = (T{x + u{x)) - R(x))V u T(x + u{x)) (3) 
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which is non-linear. Different regularisation terms S(u(x)) will lead to different non-parametric 
image registration models. Below we list three popular models selected for tests and comparisons. 
Model LC. The first one is the linear curvature model by Ei in cm ia, where 


S FMC (u) = 


(Ami) 2 + (Am 2 ) 5 


dfl. 


(4) 


This term is an approximation of the surface curvature t(ui) through the mapping ( x,y ) —>■ 
(x,y,ui(x,y)), 1= 1,2 where 


l { ui ) = V • 


Vm; 

VlVuiP + l 


Am; 


(5) 


when | Vm;| ~ 0. The Euler Lagrange equation for 0 with jjFMC ag ^g re gularisation term is 
given by a fourth order PDE 

"fA 2 u + f(u) = 0 (6) 

with boundary conditions Am; = 0, VAm; • n = 0, l = 1,2 and n the unit outward normal vector. 
The model consists of the second order derivative information of the displacement field which 
results in smoother deformations compared to those obtained using first order models based on 
elastic and diffusion energies. It is refined in m\mm with nonlinear boundary conditions. The 
affine linear transformation belongs to the kernel iS^MC( m ) j s no t the case in elastic or 

diffusion registration. 

Model MC. Next is the mean curvature model Eiia 
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MC 


(«) = / [l 

J n L 


(fc(zil)) +k(t(M 2 )) 




where k(s) = ^s 2 and t is as defined in 0. The Euler Lagrange equation for 0 with <S^4C ag 
the regularisation term is given by: 



1 

VIVmzP + I 


Vk'(t(M;)) 


Vm i • Vk'(t(M;)) 

(v^TF+1) 3 


Vm;) +/i(«) = 0, 


1 = 1,2 


(7) 


with boundary condition Vm; • n = Vi(m;) • n = 0, l = 1,2. One can use the multigrid method 
to solve equation 0 as in [3]; refer also to [2] for multi-modality image registration work. 

Model D. Finally Thirion EI; introduced the so-called demon registration method where 
every pixel in the image acts as the demons that force a pulling and pushing action in a similar 
way to what Maxwell did for solving the Gibbs paradox in thermodynamics. The original demon 
registration model is a special case of diffusion registration but it has been much studied and 
improved since 1998; see mi ns 00 nu. The energy functional for the basic demon method is 
given by 

2 

S(u) = || R(x) - T(x + u + n.) || 2 + ||'u|| 2 (8) 

where u is the current displacement field, er 2 and cr 2 account for noise on the image intensity 
and the spatial uncertainty respectively. Equation 0 can be linearised using first order Taylor 
expansion, 

2 

J(u) = ||i?(a:) - T(x + u) + Ju\\ 2 + ^-||m|| 2 (9) 

o' 2 
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where J is given by 


J= - 


Vi? + VT(a; + u) 
2 


for an efficient second order minimisation. The first order condition of © leads to the new update 
for u 


u = — 


R(x) — T(x + u ) 

lull 2 + 4 


j. 


The additional use of v for ip = exp(w) helps to achieve a nearly diffeormorphic transformation 
(mapping), where v is the stationary velocity field of the displacement field it; see [23]. It should 
be remarked that the three main steps of the model cannot be combined into a single energy 
functional. 

In a discrete setting, since the image domain S3 is a square, all variational models are discretised 
by finite differences on a uniform grid. Refer to mm- The vertex grid is defined by 


fl h = {x,„j = (Xi , yj) | 0 < i < IVi - 1 , 0 < j < N 2 - 1} 


where we shall re-use the notation T and R for discrete images of size x N%. 


3 Mathematical Background of the Gaussian curvature 

In differential geometry, the Gaussian curvature problem seeks to identify a hypersurface of R d+1 
as a graph z = u{x) over x £ C so that, at each point of the surface, the Gaussian curvature 
is prescribed. Let k(x) denote the Gaussian curvature which is a real valued function in C K d . 
The problem is modelled by the following equation 

det(.D 2 u) - k(x){1 + \Du\ 2 ) {d+2)/2 = 0 (10) 


where D is the first order derivative operator. Equation m is one of the Monge-Ampere equa¬ 
tions. For d = 2, we have 


i(x) = 


QC _ 'U'xx'Uyy 'U'xy'Uyx 

c {1 + ul + ulY 


( 11 ) 


In g], the authors define a regularisation term using the Gaussian curvature of a closed surface 
based on the Gauss-Bonnet theorem. 

Theorem 3.1 Gauss-Bonnet Theorem. For a compact C 2 surface dT,, we have 

/ K GC do = 27TY 

Jdn 

where do is the length element to the surface and x is the Euler characteristic of the surface. 

Using this Theorem, it was shown in [3] that the complete analogy of the total variation regu¬ 
larisation for surface fairing is the energy functional of the Gaussian curvature. The analogous 
term S, to the total variation of a function, that appears in the ROF model m, is given by 


S= \n{x)\do 

JdY. 


where do is the length element to the surface dT,. 
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Gaussian curvature is one of the fundamental second order geometric properties of a surface. 
According to the Gauss’s Theorema Egregium, Gaussian curvature is intrinsic. For a local iso¬ 
metric mapping / : dY —> dY' between two surfaces, Gaussian curvature remains invariant i.e. if 
p £ dY and p' £ dY', then n GC (p) = n GC (p') and the mapping / is smooth and diffeomorphic. 

We can also use a level set function to define the Gaussian curvature. Denote by 4> the zero level 
set of the surface generated through the mapping (x,y) > (x,y,u(x,y)). Then tj> = u(x,y) — z 

and Y7(j> = (u x ,u y , — 1) T where u x = and u v = The Gaussian curvature of the level set is 
given by 

= vwy (12) 

|V 0 |* ' ' 

where V0 = <j> y , (f) z ) T , |V^| = yj(j) x + (j>y + <j%, H(tj>) is the Hessian matrix and H*((j>) is the 

adjoint matrix for H(<j>). We have 


Then, 



U'XX 

'U'xy 

0 ' 


' 0 

0 

0 

H(4>) = 

'Uyx 

Uyy 

0 


0 

0 

0 


0 

0 

0 _ 


_ 0 

0 

'Uxx'U'yy ^yx^xy 


k gc 


r U / yx'U j xy 'U'xx'Uyy 
(ul+ul + l ) 2 


This is why we set k gg = —k(x) in equation (HU). We shall use \n GC \ to measure the Gaus¬ 
sian curvature as in [3] for a monotonically increasing function (since the functional should be 
nonnegative). 


4 Image Registration based on Gaussian Curvature 

Before introducing our new image registration model, we first illustrate some facts to support the 
use of Gaussian curvature. 


4.1 Advantages of a Gaussian curvature 


The total variation and Gaussian curvature. We use the volume based analysis introduced 
in m to compare two denoising models, respectively based on Gaussian curvature and the total 
variation: 




u xy u xx u yy \ 

(ul+ul + l)* ) V ’ 


(13) 


du Vu 

dt |Vu 


(14) 


Consider, for each a > 0, the time change of the volume Vt, a = {(x,y, z) \ 0 < z < \u(x, y, t) — a\} 
which is enclosed by the surface z = u(x,y,t) and the plane z = a. Assume | u(x,y,t) — a\ = 
(■ u(x,y,t ) — a)s with s either positive (s = 1) or negative (s = —1) at all points. Denote by 
Ct, a the closed curve defined by the level set u(x,y,t) = a and accordingly by dt, a the 2D region 
enclosed by Ct, a - The volume change in Vt, a in time is given by 


dr dr rHx,y,t)- a \ g , 

V = di dzdA = m I dzdA = Ft /*,. '“ (I - »■ *) - “I" 
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where dA is the area element. We now consider how V changes from evolving (fl3l) or (fill) . 

If u is the solution of equation (lldl) . then from Gauss’ theorem 

V=|/ \u(x,y,t) - a\dA = s j ^dA = s 

where da is the length element and n is the unit normal vector to the curve Ct, a which is repre¬ 
sented as n = Sj^y- Then 

,, 2 f ^ u Vu [ 

^ = S • T =-|da = / da = \c Ua \ 

Jc t , a |Vu| |Vu| J cta 

where |ct jC( | is the length of the curve c^ a . Furthermore, the volume variation in time is 


_ Vu , /' Vu 

V • -r——rdA = S 

l Vu l Jc ta |Vu| 


• nda 


\u(x, y,t + St) — a\dA 


^ t-\-St,cx 



a | dA + sSt 




a\ + St 


| C-t,a | 
|dqa | - 


dA 


where |dt, a | denotes the area of the region dt, a ■ We can see that the change in u from t, to t + St is 
proportional to the ratio . So, when this ratio is large (indicating possibly a noise presence), 
the total variation model reduces it and hence removes noise. However, important features of u 
which have a large level set ratio are removed also and so not preserved by the total variation 
model (trill . 

Using similar calculations to before for the Gaussian curvature scheme m we have 


V = dtj \u(x,y,t)-a\dA = s 


= S 


H 


i) ! 


V H 

''j Vuj ■ nda = 


*xx **yy 


{< 


i) 2 


Vu) dA 


("( 




l) 2 


^ |Vu|dcr. 


(15) 


From here, we observe that the quantity V for the subdomain v t . a is dependent on the product of 
the variation and the Gaussian curvature on the level curve. The function n in (1151) controls and 
scales the speed of the volume change in contrast to the total variation scheme where V depends 
only on the variation of the level curve. Consider a point p = (xo,yo,a) where a = u(xo,yo). 
Gaussian curvature n = k\K 2 based on two principal curvatures K\ and K 2 where Ki is the 
curvature of the level curve passing the point p and k ,2 is the curvature of the path which passes 
the point and K 2 is orthogonal to the level curve. If the Gaussian curvature on one level curve is 
zero then there is no change in V regardless of variation on the level curves. In contrast, with the 
total variation, if there is a variation in the level curve, then there is a change in V. Based on 
this observation, we believe that the Gaussian curvature model is better than the total variation 
model for preserving features on surfaces. 

The mean curvature and Gaussian curvature. The mean curvature (MC) l = (ki+« 2)/2 
is also widely used. Next, we show that, though closely related, Gaussian curvature (GC) is better 
than mean curvature for surfaces in three ways. 

First, Gaussian curvature is invariant under rigid and isometric transformations. In contrast, 
mean curvature is invariant under rigid transformations but not under isometric transforma¬ 
tions. Rigid transformations preserve distance between two points while isometric transforma¬ 
tions preserve length along surfaces and preserve angles between curves on surfaces. To illustrate 
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invariance, consider a surface 

zi(x,y) = ax 2 + by 2 , 

whose Gaussian curvature and mean curvature are respectively 

0 — (2a) (2b) (1 + 4b 2 y 2 ){2a) + (1 + Aa 2 x 2 )(2b) 

(1 + 4a 2 x 2 + 4 b 2 y 2 ) 2 ’ (1 + 4a 2 x 2 + ab 2 y 2 ) 3 / 2 

If we flip the surface upside down (isometric transformation) where z[(x 7 y) = —ax 2 — by , we 
will have the same value for the Gaussian curvature and a different value for the mean curvature. 
Thus, Gaussian curvature is invariant under isometric transformation. 

Second, Gaussian curvature can be used to localise the tip of a surface better than mean 
curvature. Consider 

Z 2 (x,y) = ~(x 2 +y 2 ). 

When we compute the mean and Gaussian curvature for the surface, we see that the value given 
by the Gaussian curvature is sharper than that of the mean curvature. The highest point of the 
Gaussian curvature is better distinguished from its neighbourhood compared to the highest point 
of the mean curvature. 

Third, Gaussian curvature can locate saddle points better than mean curvature. Take 

Z3{x,y) = -^{x 2 -y 2 ) 

as one example. We can again compute its mean and Gaussian curvatures analytically. The 
mean curvature for this surface appears complex where the largest value is not at the saddle 
point and the saddle point cannot be easily located. However, Gaussian curvature gets its highest 
value at the saddle point and is therefore able to accurately identify the saddle point within its 
neighbourhood. 

In addition to these three examples and observations, a very important fact pointed out in |4] 
states that the mean curvature of the surface is not a suitable choice for surface fairing because 
the model is not effective for preserving important features such as creases and corners on the 
surface (although the model is effective for removing noise). This is true when we refer to surface 
fairing (surface denoising) but not necessarily true for 2D image denoising. From the recent works 
done in image denoising mm, we observe several advantages of Gaussian curvature over total 
variation and mean curvature. Therefore, we conjecture that Gaussian curvature may outperform 
existing models in image registration. To our knowledge there exists no previous work on this 
topic. 

4.2 The proposed registration model 

Now we return to the problem of how to align or register two image functions T(x),R(x). Let 
the desired and unknown displacement fields between T and R be the surface map (x, y) :—>■ 
(x,y,ui(x,y)) where l = 1,2 and with u = ( ui,U 2 )• We propose our Gaussian curvature based 
image registration model as 

min J 7 (u(x)) = ]- f {T(x + u)~ R{x)) 2 dfl + ~/S GG {u(x)) (16) 

ueC 2 (Cl) Z 




where 


S GC (u(x)) = ^2s GC (ui), S GG (ui)= f 
, , Jn 


'U'l,Xy'U'l,yx 'U'l^xx'U'l^y 


( u lx + U lv 


l ) 2 


dfl. 


The above model © leads to two Euler Lagrange equations: 


7 V • 
7 V • 


4|'U2 ) a;y'a2,yx 'U j 2,xx'U j 2,yy\ 

Ml 


Vui 

Vil 2 


+ 7 V • Bi t ± + 7 V • B 12 + fi — 0 
+ 7 V • B- 2,1 + 7 ^ ■ B- 2,2 + /2 = 0 


( 17 ) 


where 


A /”l ~ U l,x + u f,y + 1) -®U — ( ( ~ 


SlUl, 


Mi 


SlUl^x 

~M7 


B 12 = 


( Sm , 


yx 


\ Mi 


Sl'lll 1 xx\ 1 c* • / \ 

J J 1 — Signal ^l,xy^l,yx ,xx^l ,yy) 


f = (/ 1 , / 2 ) T = (T{x + u) - R(x))V u T(x + «),/= 1,2 


and boundary conditions Viq ■ n = 0, Z = 1,2, where n is the normal vector at the boundary 917. 
Derivation of the resulting Euler-Lagrange equations for this model can be found in Appendix A. 
The equations are fourth order and nonlinear with anisotropic diffusion. Gradient descent type 
methods are not feasible and one way to solve them is by a geometric multigrid as in [3]- Here, 
instead, we present a fast and efficient way to solve the model (fl6l) on a unilevel grid. 


4.3 Augmented Lagrangian Method 

The augmented Lagrangian method (ALM) is often used for solving constraint minimisation 
problems by replacing the original problem with an unconstrained problem. The method is 
similar to the penalty method where the constraints are incorporated in the objective functional 
and the problem is solved using alternating minimisation of each of the sub-problems. However, 
in ALM, there are additional terms in the objective functional known as Lagrange multiplier 
terms arising when incorporating the constraints. Similar works on the augmented Lagrangian 
method in image restoration can be found in [22 E2 • 

To proceed, we introduce two new dual variables q\ and <72 where q\ = Vui(x) and <72 = 
X7u2(x). Consequently we obtain a system of second order PDEs which are more amendable to 
effective solution. 

We obtain the following refined model for Gaussian curvature image registration 
min J(u!,U 2 ,qi,q 2 ) = V(T, R,u(x)) + jS GC (q^ + jS GC (q 2 ) 

^ 1 , 142 ,< 71,<72 

s.t = Viti(£c), q 2 = Vu 2 (®) 

and further reformulate J(u\,U 2 , <Zi, ^ 72 ) to get the augmented Lagrangian functional 

C gc (u 1 ,u 2 , <7i, < 72 ; Mi, M 2 ) =^|| T(x + u{x)) - R{x)\\\ + 7<S GC (<7i) + 7 S GC (q 2 ) 

+ (£ti,< 7 i - Viii) + (7x2,92 - Vu 2 ) ( 18 ) 

+ 7 jll<?l - VU1H2 + ^||<?2 - VU2H2 

where Hh M 2 are the Lagrange multipliers, the inner products are defined via the usual integration 
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in Q and r is a positive constant. We use an alternating minimisation procedure to find the optimal 
values of ui,u 2 ,qi, q 2 and fi 1 , p 2 where the process involves only two main steps. 

Step 1. For the first step we need to update qi, q 2 for any given iti, U 2 , Mi, P 2 - The objective 
functional is given by 

mm 7 <S GC (qi) +'yS GC (q 2 ) + (mi,<?i) + (M 2 , < 22 } + t:||< 2 i - Vui || 2 + 7 -\\q 2 - Vu 2 1| 2 . 

<?1,<Z2 Z Z 

This sub-problem can be solved using the following Euler Lagrange equations: 


f f( — Ql,l)y\ , f—{qi,l)x\ 45iF»igi,2 , , r \ \ n 

-p2- ) x + (-p2- ) y ) - 7 -p3-b Ml,2 + r(qi t 2 ~ («l)y) — 0 




4<Si-Di<7i,2 


+ C—rf—“ 7 —ff— + + r(9l>1 “ (uiW “ 0 


'-(91,2) 


4 Si £ 191,1 

1.3 


(19) 


where £1 = det(Vqi) = {qi,i) x (qi, 2 )y-(qi,i)y{qi, 2 )x,^i = and 5,1 = si S n ( 

We have a closed form solution for this step, if solving alternatingly, where 


D 1 


(l|Vui|p+i ) 2 T 


9i,i = 


9 l ,2 = 


r 3(_ 7 ((^) i + (^ ¥ k)J +Mu _ r( „ I ), 

-rrf + 7 45 i^i 

r ?( - y{(^) x + (^Tr 1 ) J + - r(«l),) 


-rr? + 7 4Si£i 


Similarly, we solve 92 , 1 , 92,2 from 


f ,Y/(-92,l)y\ 

f-(q2,i)x\ 

1 7 U T\ ), + 

V r 2 A 

) (({q2a)y\ , ( 

-( 92 , 2 ) 2 , \ \ 

{ 7 U r% )„ + t 

F % )y) 


) J - 7 4 >52 ^f j2,1 + M2,1 + r(92,i - (u 2 ) x ) = 0 


( 20 ) 


, 

(WHTF+TF)' 


where £ 2 = det(Vq 2 ) = (92,i)x(92,2)j,-(92,i)j,(92,2)a, T 2 = ^-+ul x +ul y and S 2 = sign( 

Step 2. For the second step we need to update u\,u 2 for any given qi,q 2 and Mi, M 2 with 
the following functional 

min -\\T(x + u) - R{x)\\l - (mi, Vui) - (M 2 , Vu 2 ) + L<?i - Vui|| 2 + Lq 2 - V-u 2 || 2 - 

ui,u 2 Z Z Z 

Thus, we have the following Euler Lagrange equations: 


- rAui + /1 + V • Mi + rV • qi = 0 

- rAu 2 + f 2 + V • M 2 + rV ■ q 2 = 0 


( 21 ) 


with Neumann boundary conditions Viq -n = 0, l = 1,2. To solve equation (ELD, first, wc linearise 
the force term / using Taylor expansion 

fi ( u [ k+1) , u 2 k+1) ) = fl ( u i fe) , u 2 fc) ) + fl ( u [ k) , u 2 fc) )^ M / fc) + /1 ( u i fc) , u 2 k) ) Su 2 k) + ■■■ 

» fi(u[ k \u ( 2 k) ) + al k ^Su[ k) + ai^Su^ 
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where 


fi{u[ k \u^),ai t 2 = d U 2 fi(u[ k \u < ' 2 k ' l ),Su[ k) = u [ k+1 ' 1 -u[ k \ 6 u^ k) = u < ' 2 +1 ' 1 - . 

(k) (k) 

Second, we approximate a\ -[ and oy 2 with 

<j[ k ^ = [duiT{x + (du\T{x + u^)^j 

cr \ k 2 = (duiT{x + [du 2 T{x + . 

The discrete version of equation m is as follows 

N h {u h ^ k) )u h ^ k+1) = B h {u h ^ k) ) (22) 


where 


r —rC 


'n 


(u fe ’( fe )) 


' 12 


(u h 'W) 


a^(u h ^) -rC + a^ 2 (u h ^) 


N h \u <*>) = 

B h (uW) = 

L -^2 + /2W1 ' > u 2 ' 1 + V'2lV U '"' 

C is the discrete version of the Laplace operator A and Gf is the discrete version of 


r -G\ + f^(u {k \u (k) ) + a^{u^)u k ’ {k) + a^ 2 {u h 'W)u 2 ,( - k) 
~G 2 + f^\4 k) ) + a k M k) )u/ k) + a k 2 (u h ^)u h / k) 


V • fii + rV • qi, 1= 1,2. 

Third, we solve the system of equation (12^1) using a weighted pointwise Gauss Siedel method 

u K(k+i) = (1 -u)u h ’( k) +oj(N h {u {k) )y 1 B h {u ik '>) 

where u G (0, 2) and we choose w = 0.9725. 

The iterative algorithm to solve m is now summarised as follows. 

Algorithm 1 Augmented Lagrangian method for the Gaussian Curvature Image Registration. 

1. Initialise Hi = ^2 = 0 ,u(x) = 0, 7 , r\ 

2. For k = 0,1,..., IMAX 

(a) Step 1: Solve (1 1911201) for (q[ k+1 \ q^ 1 ^) with ( 111 , 112 ) = (u[ k \ ). 

(b) Step 2 : Solve (12T1) for (u^ k+1 \u' 2 +1 ' > ) with (qi,q 2 ) = (q[ k+1 \ q 2 k+1 ' > )- 

(c) Step 3: Update Lagrange multipliers. 

H { i +1) = Hi ] + r{q[ k+1) - Vw^ +1) ), /x^ fe+1) = /4*° + r(qf£ fc+1) - Vix^ +1) ) 

3. End for. 


5 Numerical Results 

We use two numerical experiments to examine the efficiency and robustness of the Algorithm 
Q] on a variety of deformations. To judge the quality of the alignment we calculate the relative 
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reduction of the similarity measure 


T>(T, R, u) 
V{T, R) 


and the minimum value of the determinant of the Jacobian matrix J of the transformation, 
denoted as T 


J = 


1 + Ml,, 

U2,x 


u l,y 

1 + U 2 ■ 


JF = min (det(J)). 


(23) 


We can observe that when T > 0, the deformed grid is free from folding and cracking. 

All experiments were run on a single level. Experimentally, we found that r £ [0.02,2] works 
well for several types of images. As for the stopping criterion, we use tol = 0.001 for the residual of 
the Euler-Lagrange equations m-m and the maximum number of iterations is 30. Experiments 
were carried out using Matlab R2014b with Intel(R) core (TM) i7-2600 processor and 16G RAM. 


5.1 Test 1: A Pair of Smooth X-ray Images 

Images for Test 1 are taken from m where X-ray images of two hands of different individuals 
need to be aligned. The size of images are 128 x 128 and the recovered transformation is expected 
to be smooth. The scaled version of the transformation and the transformed template image is 
given in Figures [1] (d) and (e) respectively. The transformation is smooth and the model is able 
to solve such a problem. For comparison, the transformed template images for the diffeormorphic 



(a )T 



(d) x + u(x) 


(b) R 



(e) T(x + u(x)) 


(c) T - R 



(f) T(x + u(x)) — R 


Figure 1: Test 1 (X-ray hand). Illustration of the effectiveness of Gaussian curvature with smooth 
problems. On the top row, from left to right: (a) template, (b) reference and (c) the difference 
before registration. On the bottom row, from left to right: (d) the transformation applied to a 
regular grid, (e) the transformed template image and (f) the difference after registration. As can 
be seen from the result (e) and the small difference after registration (f), Gaussian curvature is 
able to solve smooth problems. 


demon method, linear, mean and Gaussian curvatures are shown in Figures [2] (a), (b), (c) and 
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(d) respectively. We can observe that there are some differences of these images inside the red 
boxes where only Gaussian curvature delivering the best result of the features inside the boxes. 



(a) Model D 


(b) Model LC 



(c) Model MC (d) Gaussian curvature 

Figure 2: Test 1 (X-ray hand). Comparison of Gaussian curvature with competing methods. The 
transformed template image using (a) Model D, (b) Model LC, (c) Model MC and (d) Gaussian 
curvature. Note the difference of these three images inside the red boxes. 


Measure 

Model D 

Model LC 

Model MC 

GC 


ML 

SL 

ML 

SL 

SL 

SL 

7 

1.6 

1.6 

0.1 

0.5 

0.0001 

0.0001 

Time (s) 

15.19 

186.48 

84.33 

12.98 

275.3 

953.15 

£ 

0.1389 

0.1229 

0.0720 

0.3780 

0.0964 

0.0582 

T 

0.0600 

0.1082 

0.3894 

0.1973 

0.6390 

0.3264 


Table 1: Quantitative measurements for all models for Test 1. ML and SL stand for multi and 
single level respectively. 7 is chosen as small as possible such that T > 0 for all methods. T > 0 
indicates the deformation consists of no folding and cracking of the deformed grid. We can see 
that the smallest value of e is given by Gaussian curvature (GC). 


We summarised the results for Test 1 in Table [1] where ML and SL stand for multi and single 
level respectively. For all models, 7 is chosen as small as possible such that T > 0. We can see 
that the fastest model is the diffeormorphic demon, followed by linear and mean curvature. The 
current implementation for Gaussian curvature is on single level and the model uses augmented 
Lagrangian method which has four dual variables and four lagrange multipliers terms. Thus, it 
requires more computational time than the other models. 


13 


















5.2 Test 2: A Pair of Brain MR Images 


We take as Test 2 a pair of medical images of size 256 x 256 from the Internet Brain Segmentation 
Repository (IBSR) http: //www. cma. mgh. harvard. edu/ibsr where 20 normal MR brain images 
and their manual segmentations are provided. We choose the particular pair of individuals with 
different sizes of ventricle to illustrate a large deformation problem. Figure [3] shows the test images 
and the registration results using Gaussian curvature model. We can see that the model is able 
to solve real medical problems involving large deformations, which is particularly important for 
atlas construction in medical applications. Figure U shows the transformed template images for 



(d) x + u{x) 


(e) T(x + u(x)) 



(c )T- R 



(f) T(x + u{x)) -R 


Figure 3: Test 2: A pair of Brain MR images. Illustration of the effectiveness of Gaussian cur¬ 
vature with real medical images. On the top row, from left to right: (a) template, (b) reference 
and (c) the difference before registration. On the bottom row, from left to right: (d) the trans¬ 
formation applied to a regular grid, (e) the transformed template image and (f) the difference 
after registration. As can be seen from the result (e) and the small difference after registration 
(f), Gaussian curvature can be applied to real medical images and is able to obtain good results. 


all four methods. We can see that Gaussian curvature gives the best result inside the red boxes 
in comparison with the diffeomorphic demon, the linear and mean curvature models as depicted 
in Figure [4] (d). 

The values of the quantitative measurements for Test 2 are recorded in Table [2] where the 
lowest values of e are given by the Gaussian curvature model indicating higher similarity between 
the transformed template result and the reference image. However, our proposed model required 
more time than the other models since the model consists more variables than the others. 

A brief summary. The linear curvature model is relatively easy to solve, based on ap¬ 
proximation of the mean curvature. The mean curvature model for image registration is highly 
nonlinear, making it challenging to solve. The Gaussian curvature resembles the mean curvature 
in many ways, though different, but its model appears to deliver better quality than the mean 
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(a) Model D 


(b) Model LC 




(c) Model MC (d) Gaussian curvature 

Figure 4: Test 2: A pair of Brain MR images. Comparison of Gaussian curvature with competing 
methods. The transformed template image using (a) Model D, (b) Model LC, (c) Model MC, 
and (d) Gaussian curvature. Notice the differences of these three images inside the red boxes. 
Considerably more accurate results are obtained, particularly within these significant regions, by 
employment of the Gaussian curvature model. 


Measure 

Model D 

Model LC 

Model MC 

GC 


ML 

SL 

ML 

SL 

SL 

SL 

7 

1.2 

1.4 

0.16 

2.0 

0.0001 

0.0001 

Time (s) 

23.89 

209.00 

275.04 

35.70 

830.22 

1053.7 

£ 

0.2004 

0.7580 

0.1128 

0.4283 

0.1998 

0.1062 

T 

0.0277 

0.0387 

0.3157 

0.0148 

0.8240 

0.0138 


Table 2: Quantitative measurements for all models for Test 2. ML and SL stand for multi and 
single level respectively. 7 is chosen to be as small as possible such that T > 0 for all models. 
T > 0 indicates the deformation consists of no folding and cracking of the deformed grid. We 
can see that the smallest value of e is given by Gaussian curvature (GC). 


curvature. The diffeormorphic demon model is equivalent to the second order gradient descent 
on the SSD as shown in m The model is only limited to mono-modality images and it is not 
yet applicable to multi-modality images. Our Gaussian curvature model however can be easily 
modified to work with multi-modality images by replacing the SSD by a mutual information or 
normalised gradient fields based regularizer; an optimal solver is yet to be developed. We show 
one example of extension to a pair of multi-modality images in Figj5] 


6 Conclusions 

We have introduced a novel regularisation term for non-parametric image registration based on 
the Gaussian curvature of the surface induced by the displacement held. The model can be 
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(c) T (x + u(x)) (d) x + u(x) 

Figure 5: Results of Gaussian curvature image registration for multi-modality images. The model 
is able to register multi-modality images with mutual information as the distance measure. 

effectively solved using the augmented Lagrangian and alternating minimisation methods. For 
comparison, we used three models: the linear curvature [Bj, the mean curvature [3] and the demon 
algorithm [51] for mono-modality images. Numerical experiments show that the proposed model 
delivers better results than the competing models. 


Appendix A — Derivation of the Euler-Lagrange Equations 

Let q i = u x and q 2 = u y \ then we can write the Gaussian curvature regularisation term as 


S GC ( qi ,q 2 )= f 
J n 


Q\,xQ2,y ~ qi,yq2,x 


{l + ql + ql) 2 


dxdy. 


From the optimality condition 


and ^£ GC 0?i;92 +£ 2 ^ 2 ) 


S<n -^-= 0, then ^-S GC (q 1 +e 1 i Pl ,q 2 ) 

= 0. In details, 


= 0 


ei=0 


€ 2=0 


— [ 

de 1 J n 


(gi + eiipi) x q 2lV - (qi + tx<Px) v q2,x 
(1 + (qx +e 1 yi) 2 +g|) 2 


dxdy 


e —0 


[ s d 

'{qx + ei<Pi)xQ2,y - (qi + ei<Pi) v q 2 ,x' 



Index 

(1 + {qi + eiv?i) 2 + ql) 2 

UitAsUi y 

6=0 


(24) 
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where S = sign a 2 ’* ). From 


in 


Q P l,xQ2,y ( Pl,yQ2,x , \l A (A , 2 , 2\-3\ 

5 (1 + q 2 + g 2)2 + (gMg 2 ,y ~ gi, y g2,x)(~4^lgl(l + q 1 + ? 2 ) ) 


dxdy 


Sifii,xQ2,y Sipi, y q 2 , x ASDqiipi 


Jn F 2 r 2 T3 

where T = 1 + q\ + q\, D = qi, x q 2 , y - qi, v q 2 ,x 


dxdy = 0 , 


Using the Green theorem f 9Q <t>u> ■ nds — f n (f>div(u>)dxdy = f n V<j> ■ udxdy, we have, 


/' S(p-L, x q 2 ,y Stfii y q 2 x , , /' (Sq 2 .y Sq 2 x \ /' f Sq 2 , y Sq 2; 

L — -f 1 ildv = L Vl {~ r^- — j - - J a 


= 0 


where <f> = yq, lS p 2 ,a; ^. Setting the boundary integral to zero, then we derive 


/' j. CSg 2 , y Sq 2 , x \ 

J=°- 

Finally, we use the fundamental lemma of calculus of variation to get: 

„ fSq 2 , y Sq 2 , x \ ASDqi = 

V \ F 2 ’ T 2 / F 3 


Similarly, for ^S GC (q 1 ,q 2 + 62 ^ 2 ) 


0 , we finally obtain equation (ED. 


□ 
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